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Abstract 

In the earlier work on the development of a model-independent data analysis method for 
determining the mass of Weakly Interacting Massive Particles (WIMPs) by using measured 
recoil energies from direct Dark Matter detection experiments directly, it was assumed that 
the analyzed data sets are background-free, i.e., all events are WIMP signals. In this article, 
as a more realistic study, we take into account a fraction of possible residue background 
events, which pass all discrimination criteria and then mix with other real WIMP-induced 
events in our data sets. Our simulations show that, for the determination of the WIMP 
mass, the maximal acceptable fraction of residue background events in the analyzed data 
sets of O(50) total events is ~ 20%, for background windows of the entire experimental 
possible energy ranges, or in low energy ranges; while, for background windows in relatively 
higher energy ranges, this maximal acceptable fraction of residue background events can not 
be larger than ~ 10%. For a WIMP mass of 100 GeV with 20% background events in the 
windows of the entire experimental possible energy ranges, the reconstructed WIMP mass 
and the la statistical uncertainty are ~ 97 GeV ^35% (~ 94 GeV ^33% for background-free 
data sets). 



1 Introduction 



Currently, direct Dark Matter detection experiments searching for Weakly Interacting Massive 
Particles (WIMPs) are one of the promising methods for understanding the nature of Dark 
Matter and identifying them among new particles produced at colliders as well as reconstructing 
the (sub) structure of our Galactic halo [1, 2, 3, 4]. In order to determine the mass of halo 
WIMPs without making any assumptions about their density near the Earth or their velocity 
distribution nor knowing their scattering cross section on nucleus, a model-independent method 
by combining two experimental data sets with two different target nuclei has been developed 
[5, 6]. This method builds on the earlier work on the reconstruction of the (moments of the) 
one-dimensional velocity distribution function of halo WIMPs, fi(v), by using data from direct 
detection experiments [7]. 

In the analysis of reconstructing fi(v), the moments of the WIMP velocity distribution 
function can be determined from experimental data directly with an unique input information 
about the WIMP mass m x . Hence, one can simply require that the values of a given moment 
of fi(v) determined by two experiments agree 1 . This leads to a simple analytic expression 
for determining m x [5, 6], where each moment can in principle be used. Additionally, under 
the assumptions that the spin-independent (SI) WIMP-nucleus interaction dominates over the 
spin-dependent (SD) one and the SI WIMP coupling on protons is approximately the same as 
that on neutrons, a second analytic expression for determining m x has been derived [6]. Finally, 
by combining the first estimators for different moments with each other and with the second 
estimator, one can yield the best-fit WIMP mass as well as minimize its statistical uncertainty. 

In the work on the development of the model-independent data analysis procedure for the 
determination of the WIMP mass, it was assumed that the analyzed data sets are background- 
free, i.e., all events are WIMP signals. Active background discrimination techniques should make 
this condition possible. For example, the ratio of the ionization to recoil energy, the so-called 
"ionization yield", used in the CDMS-II experiment provides an event-by-event rejection of 
electron recoil events to be better than misidentification [8]. By combining the "phonon 
pulse timing parameter", the rejection ability of the misidentified electron recoils (most of them 
are "surface events" with sufficiently reduced ionization energies) can be improved to be < 1CT 6 
for electron recoils [8]. Moreover, as demonstrated by the CRESST collaboration [9], by means 
of inserting a scintillating foil, which causes some additional scintillation light for events induced 
by a-decay of 210 Po and thus shifts the pulse shapes of these events faster than pulses induced 
by WIMP interactions in the crystal, the pulse shape discrimination (PSD) technique can then 
easily distinguish WIMP-induced nuclear recoils from those induced by backgrounds 2 . 

However, as the most important issue in all underground experiments, the signal identification 
ability and possible residue background events which pass all discrimination criteria and then mix 
with other real WIMP-induced events in our data sets should also be considered. Therefore, 
in this article, as a more realistic study, we take into account different fractions of residue 
background events mixed in experimental data sets and want to study how well the model- 
independent method could reconstruct the input WIMP mass by using these "impure" data sets 
and how "dirty" these data sets could be to be still useful. 

The remainder of this article is organized as follows. In Sec. 2 we review the recoil spectrum 
of elastic WIMP-nucleus scattering and introduce two kinds of background spectrum used in 

1 Note that, as demonstrated and discussed in Ref. [6], this condition requires an algorithmic procedure for 
matching the maximal cut-off energies of the analyzed data sets. 

2 For more details about background discrimination techniques and status in currently running and projected 
direct detection experiments, see e.g., Refs. [10, 11, 12] 



our simulations. In Sec. 3 we first review briefly the model-independent method for the deter- 
mination of the WIMP mass. Then we show numerical results of the reconstructed WIMP mass 
by using mixed data sets with different fractions of residue background events based on Monte 
Carlo simulations. We conclude in Sec. 4. Some technical details will be given in an appendix. 



2 Signal and background spectra 

In this section we first review the recoil spectrum of elastic WIMP-nucleus scattering. Then 
we introduce two forms of background spectrum which will be used in our simulations. Some 
numerical results of the measured energy spectrum superposed by the WIMP scattering and 
background spectra will also be discussed. 

2.1 Elastic WIMP— nucleus scattering spectrum 

The basic expression for the differential event rate for elastic WIMP-nucleus scattering is given 
by [3]: 



dQ ^ ^ Iv mi 



fi(v) 
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Here R is the direct detection event rate, i.e., the number of events per unit time and unit mass 
of detector material, Q is the energy deposited in the detector, F(Q) is the elastic nuclear form 
factor, f\(v) is the one-dimensional velocity distribution function of the WIMPs impinging on 
the detector, v is the absolute value of the WIMP velocity in the laboratory frame. The constant 
coefficient A is defined as 

A=j^, (2) 

where p is the WIMP density near the Earth and a Q is the total cross section ignoring the form 
factor suppression. The reduced mass m r N is defined by 

™ r ,N = ; , (3) 

m x + mN 

where m x is the WIMP mass and m N that of the target nucleus. Finally, t> min is the minimal 
incoming velocity of incident WIMPs that can deposit the energy Q in the detector: 



Umin = ol^Q (4) 
with the transformation constant 



and i> max is the maximal WIMP velocity in the Earth's reference frame, which is related to 
the escape velocity from our Galaxy at the position of the Solar system, v esc > 600 km/s. 
Note that, as will be shown below, the Earth's velocity relative to the Galactic halo is time- 
dependent, and considering the random motion of WIMPs in the Galaxy, the relation between 
the one-dimensional cut-off v max and the three-dimensional one v esc is thus rather complicated. 
Nevertheless, it is unlike to affect the event rate as well as the results shown in this article 
significantly. In the literature, for simplicity and practical uses, t> max is often set as oo (e.g., 
[13, 14, 15]). 



2.1.1 One dimensional WIMP velocity distribution function 



The simplest semi-realistic model halo is a spherical isothermal Maxwellian halo. More realisti- 
cally, one has to take into account the orbital motion of the Solar system around the Galaxy as 
well as that of the Earth around the Sun. The one-dimensional velocity distribution function of 
this shifted Maxwellian halo can be expressed as [2, 3, 7] 

1 



/l,Bh(u) = —J= 



(6) 



-(v-v c ) 2 /v 2 _ p -(v+v c ) 2 /v 2 

111 \V C V ,' ' 

Here vq — 220 km/s is the orbital velocity of the Sun in the Galactic frame, and v e is the Earth's 
velocity in the Galactic frame [16, 3, 4]: 

'27i(t-t p y ~ 



v e (t) = v 



1.05 + 0.07 cos 



1 yr 



(7) 



Lp — 



June 2nd is the date on which the velocity of the Earth relative to the WIMP halo is 
maximal. Substituting Eq. (6) into Eq. (1), an analytic form of the integral over the velocity 
distribution function can be given as [17] 
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Here erf(x) is the error function, defined as 
2 
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On the other hand, for practical, numerical uses, an approximate form of the integral over 
fi(v) was introduced as [2] 
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where c and C\ are two fitting parameters of order unity. Not surprisingly, their values depend 
on the Galactic orbital and escape velocities, the target nucleus, the threshold energy of the 
experiment, as well as on the mass of incident WIMPs. Note that, the characteristic energy 
Qch = ciVq/cx 2 and thus the shape of the recoil spectrum depend highly on the WIMP mass: 
for light WIMPs (m x ton), Qch oc and the recoil spectrum drops sharply with increasing 



recoil energy, while for heavy WIMPs 
flatter. 



3> ton), Qch ~ const, and the spectrum becomes 



2.1.2 Spin-independent WIMP-nucleus cross section 

In most theoretical models, the spin-independent (SI) WIMP-nucleus interaction with an atomic 
mass number A > 30 dominates over the spin-dependent (SD) one [3, 4]. Additionally, for the 
lightest supersymmetric neutralino which is perhaps the best motivated WIMP candidate [3, 4], 
and for all WIMPs which interact primarily through Higgs exchange, the SI scalar coupling is 
approximately the same on both protons p and neutrons n, the "pointlike" cross section a in 
Eq. (2) can thus be written as 

*o = W^) 2 .S, do) 



where 
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is the SI WIMP-proton cross section, f p is the effective XXPP four-point coupling, A is the 
atomic mass number of the target nucleus, and m r p is the reduced mass of the WIMP mass m x 
and the proton mass m p . 

For the SI WIMP-nucleus cross section, an analytic form for the elastic nuclear form factor, 
inspired by the Woods-Saxon nuclear density profile, has been suggested by Engel as [18, 3, 4] 3 
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Here j\(x) is a spherical Bessel function, q = \j2m^Q is the transferred 3-momentum, given as 
a function of the recoil energy transferred from the incident WIMP to the target nucleus, Q, 



and the mass of the target nucleus, m N ; R 1 = \jR\ — 5s 2 is the effective nuclear radius with 
Ra — 1.2 A 1 / 3 fm and the nuclear skin thickness s ~ 1 fm. 

2.2 Background spectrum 

For our simulations with residue background events, two forms of background spectrum are 
considered. The simplest choice for the background spectrum is the constant spectrum: 

More realistically, inspired by Ref. [14], we introduce a target-dependent exponential spectrum 
given by 

Here Q is the recoil energy, A is the atomic mass number of the target nucleus. The power index 
of A, 0.6, is an empirical constant, which has been chosen so that the exponential background 
spectrum is somehow similar to, but still different from the expected recoil spectrum of the target 
nuclei; otherwise, there is in practice no difference between the WIMP scattering and background 
spectra. Note that, among different possible choices (e.g., the exponential form used in Ref. [14]), 
we use in our simulations the atomic mass number A as the simplest, unique characteristic 
parameter in the general analytic form (14) for defining the residue background spectrum for 
different target nuclei. However, it does not mean that the (superposition of the real) background 
spectra would depend simply/primarily on A or on the mass of the target nucleus, m^. In other 
words, it is practically equivalent to use expression (14) or (dR/dQ)\, giCX = e _Q / 13 - 5 keV directly 
for a 76 Ge target. 

Note also that, firstly, two forms of background spectrum given in Eqs. (13) and (14) are 
rather naive; however, since we consider here only a few residue background events induced by 
perhaps two or more different sources, pass all discrimination criteria, and then mix with other 
WIMP-induced events in our data sets of O(50) total events, exact forms of different background 
spectra are actually not very important and these two spectra, in particular, the exponential 
one, should practically not be unrealistic 4 . Secondly, for using the maximum likelihood analysis 

3 Other commonly used analytic forms for the nuclear form factor for the SI WIMP-nucleus cross section can 
be found in Ref. [17]. 

4 Other (more realistic) forms for background spectrum (perhaps also for some specified targets/experiments) 
can be tested on the AMIDAS website [19, 20]. 



to determine the WIMP mass, as described in Refs. [21, 14, 13], a prior knowledge about the 
WIMP scattering spectrum and eventually about the background spectrum is essential [14]. 
In contrast, as demonstrated in Ref. [6] and will be reviewed in the next section, the model- 
independent data analysis procedure requires only measured recoil energies (induced mostly by 
WIMPs and occasionally by background sources) from two experimental data sets with different 
target nuclei. Therefore, for applying this method to future real data from direct detection 
experiments, the prior knowledge about (different) background source(s) is not required at all. 

2.3 Measured energy spectrum 

In Figs. 1 we show measured energy spectra (solid red histograms) for a 76 Ge target with six dif- 
ferent WIMP masses: 10, 25, 50, 100, 250, and 500 GeV based on Monte Carlo simulations. The 
dotted blue curves are the elastic WIMP-nucleus scattering spectra for the shifted Maxwellian 
velocity distribution given in Eq. (6) with t> = 220 km/s, v c = 1.05 v , 5 and f esc = 700 km/s 
and the Woods-Saxon elastic nuclear form factor in Eq. (12). The dashed green curves are the 
exponential background spectra given in Eq. (14), which have been normalized so that the ratios 
of the areas under these background spectra to those under the (dotted blue) WIMP scattering 
spectra are equal to the background-signal ratio in the whole data sets (i.e., 20% backgrounds 
to 80% signals shown in Figs. 1). The experimental threshold energy has been assumed to be 
negligible and the maximal cut-off energy is set as 100 keV. 5,000 experiments with 500 total 
events on average in each experiment have been simulated. 

The measured energy spectra (solid red histograms) shown in Figs. 1 are averaged over 
the simulated experiments. Five bins with linearly increased bin widths have been used for 
binning generated signal and background events. As argued in Ref. [7], for reconstructing the 
one-dimensional WIMP velocity distribution function, this unusual, particular binning has been 
chosen in order to accumulate more events in high energy ranges and thus to reduce the statistical 
uncertainties in high velocity ranges. However, as we will show later, for the determination of 
the WIMP mass, one needs either events in the first energy bin or all events in the whole data 
set. Hence, there is in practice no difference between using an equal bin width for all bins or the 
(linearly) increased bin widths. 

Note here that, firstly, the possible energy ranges in which residue background events exist 
(the background windows) have been assumed to be the same as the entire experimental possible 
energy ranges (e.g., between and 100 keV for simulations shown in Figs. 1). Secondly, the actual 
numbers of signal and background events in each simulated experiment are Poisson-distributed 
around their expectation values independently. This means that, for example, for simulations 
shown in Figs. 1 we generate 400 (100) events on average for WIMP signals (backgrounds) 
and the total event number recorded in one experiment is then the sum of these two numbers. 
Thirdly, for the simulations demonstrated here as well as in the next section, we assumed that 
all experimental systematic uncertainties as well as the uncertainty on the measurement of the 
recoil energy could be ignored. The energy resolution of most existing detectors is so good that 
its error can be neglected compared to the statistical uncertainty for the foreseeable future with 
pretty few events. 

In Figs. 1 it can be found that, as mentioned earlier, the shape of the WIMP scattering 
spectrum depends highly on the WIMP mass: for light WIMPs (m x < 50 GeV), the recoil spectra 
drop sharply with increasing recoil energies, while for heavy WIMPs (m x > 100 GeV), the spectra 
become flatter. In contrast, the exponential background spectra shown here depend only on the 

5 Thc time dependence of the Earth's velocity in the Galactic frame, the second term of v c (t) in Eq. (7), has 
been ignored. 
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Figure 1: Measured energy spectra (solid red histograms) for a 76 Ge target with six different 
WIMP masses: 10, 25, 50, 100, 250, and 500 GeV. The dotted blue curves are the elastic WIMP- 
nucleus scattering spectra for the shifted Maxwellian velocity distribution and the Woods-Saxon 
elastic nuclear form factor; whereas the dashed green curves are the exponential background 
spectra normalized to fit to the chosen background ratio, which has been set as 20% here. 
The experimental threshold energy has been assumed to be negligible and the maximal cut-off 
energy is set as 100 keV. The background windows have been assumed to be the same as the 
experimental possible energy ranges. 5,000 experiments with 500 total events on average in each 
experiment have been simulated. See the text for further details. 
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Figure 2: As in Figs. 1, except that the constant background spectrum in Eq. (13) has been 
used. See the text for further details. 



target mass and are rather flatter (sharper) for light (heavy) WIMP masses compared to the 
WIMP scattering spectra. This means that, once input WIMPs are light (heavy), background 
events would contribute relatively more to high (low) energy ranges, and, consequently, the 
measured energy spectra would mimic scattering spectra induced by heavier (lighter) WIMPs. 

As a comparison, in Figs. 2 we generate background events with the constant spectrum 
given in Eq. (13). It can be seen clearly that, since the background spectrum now is flatter for 
all WIMP masses, background events contribute always relatively more to high energy ranges, 
and, therefore, the measured energy spectra would always mimic scattering spectra induced by 
heavier WIMPs. 



3 Reconstruction of the WIMP mass 

In this section we first review the model-independent method for determining the WIMP mass 
introduced in Refs. [5, 6]. Then we demonstrate some numerical results of the reconstructed 
WIMP mass by using mixed data sets from WIMP signals and background events based on 
Monte Carlo simulations. 

3.1 Model— independent determination of the WIMP mass 

Here we review briefly the model-independent data analysis procedure for the determination 
of the WIMP mass by using two experimental data sets with different target nuclei. Detailed 
derivations and discussions can be found in Refs. [5, 6]. 

3.1.1 Basic expressions for determining the WIMP mass 

In the earlier work [7], it was found that the normalized one-dimensional velocity distribution 
function of incident WIMPs can be solved from Eq. (1) directly and, consequently, its generalized 
moments can be estimated by [6] 
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is an estimated value of the measured recoil spectrum (dR/dQ) cxpt (before the normalization by 
the exposure S) at Q = Q m in, and I n (Qmm, Qmax) can be estimated through the sum: 
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By requiring that the values of a given moment of fi(v) estimated by Eq. (15) from two 
experiments with different target nuclei, X and Y, agree, m x appearing in the prefactor a n on 
the right-hand side of Eq. (15) has been solved as [5]: 



m 



X\(v n ) 



^/m x m Y - m x {TZ n ,x/Tln,Y) 
Kn,x /Hn,Y ~ \Jm x /m Y 



where 



n,X 



2Q^i /2 rx(g r 



1 ,x)/i 7 l(Qmin,x) + (n+l) 



2Qmin, 

X r x(Qmin,x)/F x (Q 

min,X ) + h,x 



l/n 



(19) 



and TZ n ,Y can be defined analogously. Here n^O, mrx,Y) and F/x,y)(Q) are the masses and the 
form factors of the nucleus X and Y, respectively, and rrx,Y)(Qimn,(x,Y)) refer to the counting 
rates for detectors X and Y at the respective lowest recoil energies included in the analysis. Note 
that, firstly, the general expression (18) can be used either for spin-independent or for spin- 
dependent scattering, one only needs to choose different form factors under different assumptions. 
Secondly, the form factors in the estimate of I HyX and l n y using Eq. (17) are also different. 

On the other hand, by using the theoretical prediction that the SI WIMP-nucleus cross 
section dominates, and the fact that the integral over the one-dimensional WIMP velocity 
distribution on the right-hand side of Eq. (1) is the minus-first moment of this distribution, 
which can be estimated by Eq. (15) with n = —1, one can easily find that [6] 
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Note that the exposure of the experiment, £ , appears in the denominator. Since the unknown 
factor po|/p| 2 on the left-hand side above is identical for different targets, it leads to a second 
expression for determining m x [6]: 
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and similarly for IZ^y. 



3.1.2 x 2 fitting 

In order to yield the best-fit WIMP mass as well as to minimize its statistical uncertainty by 
combining the estimators for different n in Eq. (18) with each other and with the estimator in 
Eq. (21), a x 2 function has been introduced as [6] 
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the other n max + 2 functions /^y can be defined analogously. Here n max determines the highest 
moment of fi(v) that is included in the fit. The fi are normalized such that they are dimension- 
less and very roughly of order unity in order to alleviate numerical problems associated with the 
inversion of their covariance matrix. Note that the first n max + 1 fit functions depend on m x only 
through the overall factor a and that m x in Eqs. (24a) and (24b) is now a fit parameter, which 
may differ from the true value of the WIMP mass. Finally, C in Eq. (23) is the total covariance 
matrix. Since the X and Y quantities are statistically completely independent, C can be written 
as a sum of two terms: 



&j = cov (f i<x , fj,x) + cov (f i>Y , fj, Y ) ■ 



(25) 



3.1.3 Matching the cut— off energies 

The basic requirement of the expressions for determining m x given in Eqs. (18) and (21) is that, 
from two experiments with different target nuclei, the values of a given moment of the WIMP 
velocity distribution estimated by Eq. (15) should agree. This means that the upper cuts on 
fi(v) in two data sets should be (approximately) equal 6 . Since t> cut = ay/Q max , it requires that 
[6] 

Qmax,y=( — ) Qmax,X ■ (26) 

Note that a defined in Eq. (5) is a function of the true WIMP mass. Thus this relation for 
matching optimal cut-off energies can be used only if m x is already known. One possibility to 
overcome this problem is to fix the cut-off energy of the experiment with the heavier target, 
minimize the x 2 { m x) function defined in Eq. (23), and then estimate the cut-off energy for the 
lighter nucleus by Eq. (26) algorithmically [6]. 



3.2 Reconstructing m x by using data sets with background events 

In this subsection we show some numerical results of the reconstruction of the WIMP mass 
with mixed data sets from WIMP-induced and background events by means of the model- 
independent method described in the previous subsection. The upper and lower bounds on the 
reconstructed WIMP mass are estimated from the requirement that \ 2 exceeds its minimum by 
l. 7 As in Ref. [6], 28 Si and 76 Ge have been chosen as two target nuclei. The scattering cross 
section o"o in Eq. (2) has been assumed to be dominated by the spin-independent WIMP-nucleus 
interaction. The experimental threshold energies of two experiments have been assumed to be 
negligible and the maximal cut-off energies are set the same as 100 keV. 2 x 5,000 experiments 
have been simulated. In order to avoid large contributions from very few events in high energy 
ranges to the higher moments [7], only the moments up to n max = 2 were included in the x 2 fit. 

6 Here the threshold energies have been assumed to be negligibly small. 

7 Note that, rather than the mean values, the (bounds on the) reconstructed WIMP mass are always the 
median values of the simulated results. 
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Figure 3: The reconstructed WIMP mass and the lower and upper bounds of the la statistical 
uncertainty with mixed data sets from WIMP-induced and background events as functions of 
the input WIMP mass. 28 Si and 76 Ge have been chosen as two target nuclei. The background 
ratios shown here are no background (dashed green curves), 10% (long-dotted blue curves), 20% 
(solid red curves), and 40% (dash-dotted cyan curves) background events in the whole data sets 
in the experimental energy ranges between and 100 keV. Each experiment contains 50 total 
events on average before cuts on Q max for the experiments with the Si target; all of these events 
are treated as WIMP signals. Other parameters are as in Figs. 1. See the text for further details. 

3.2.1 With the exponential background spectrum 

Fig. 3 shows the reconstructed WIMP mass and the lower and upper bounds of the la statistical 
uncertainty with mixed data sets from WIMP-induced and background events as functions of 
the input WIMP mass. As in Figs. 1, the exponential background spectrum has been used 
and the background windows are set as the same as the experimental possible energy ranges, 
i.e., between and 100 keV for both experiments. The background ratios shown here are no 
background (dashed green curves), 10% (long-dotted blue curves), 20% (solid red curves), and 
40% (dash-dotted cyan curves) background events in the whole data sets. Each experiment 
contains 50 total events on average before cuts on Q max for the experiments with the Si target. 
Remind that all events recorded in our data sets are treated as WIMP signals in the analysis, 
although statistically we know that a fraction of these events could be backgrounds. 

It can be seen clearly that, for light WIMP masses (m x < 100 GeV), the larger the fraction 
of background events in the data sets, the heavier the reconstructed WIMP masses as well as 
the statistical uncertainty intervals. This is caused directly by the background contribution 
to high energy ranges shown in Figs. 1. As discussed in Sec. 2.3, the background spectrum 
is relatively flatter compared to the scattering spectrum induced by light WIMPs, and the 
energy spectrum of all recorded events would thus mimic a scattering spectrum induced by 



WIMPs with a relatively heavier mass. Not surprisingly, the larger the background ratio, the 
more the background contribution to high energy ranges, and, consequently, the more strongly 
overestimated the reconstructed WIMP masses as well as the statistical uncertainty intervals. 

In contrast, for heavy WIMP masses (m x > 100 GeV), Fig. 3 does not show very clearly but a 
tendency 8 that the larger the fraction of background events, the lighter the reconstructed WIMP 
masses as well as the statistical uncertainty intervals. This is now caused by the background 
contribution to low energy ranges shown in Figs. 1. As discussed in the previous section, the 
background spectrum is relatively sharper compared to the scattering spectrum induced by heavy 
WIMPs, and the energy spectrum of all recorded events would thus mimic a scattering spectrum 
induced by WIMPs with a relatively lighter mass. Moreover, the larger the background ratio, the 
more the background contribution to low energy ranges, and, consequently, the more strongly 
underestimated the reconstructed WIMP masses as well as the statistical uncertainty intervals. 

Nevertheless, from Fig. 3 it can be found that, with ~ 20% residue background events in 
the analyzed data sets, the true values of the WIMP mass can still fall in the middle of the 
la statistical uncertainty band and one could thus in principle reconstruct the WIMP mass 
pretty well; if WIMPs are light (m x < 200 GeV), the maximal acceptable fraction of residue 
background events could even be as large as ~ 40%. For a WIMP mass of 100 GeV with 20% 
background events in the data sets, the reconstructed WIMP mass and the statistical uncertainty 
are ~ 97 GeV ^35%, compared to ~ 94 GeV ~^H% for background-free data sets; for a lighter 
WIMP mass of 50 GeV, the reconstructed WIMP mass and the statistical uncertainty change 
from ~ 48 GeV+^g (background-free), to ~ 54 GeV+^ (20% background), and ~ 61 GeV^S 
(40% background). 

On the other hand, considering different efficiencies of discrimination ability against different 
background sources in different energy ranges in different experiments, in Figs. 4 we shrink the 
background window in each experiment to a relatively lower range between and 50 keV (upper) 
and a relatively higher range between 50 and 100 keV (lower) 9 . Since our background spectrum 
is exponential, for the case shown in Fig. 3, only very few background events could be observed in 
the energy range between 50 and 100 keV. Hence, for the case with the background window only 
in the low energy range, not surprisingly, the results of the reconstructed WIMP mass shown in 
the upper frame of Figs. 4 should not differ very much from those shown in Fig. 3. However, 
due to the little bit more contribution to the low energy range from background events, all the 
reconstructed WIMP masses shown here are somehow lighter than those shown in Fig. 3. Hence, 
with ~ 20% residue background events in low experimental possible energy ranges, one could in 
principle reconstruct the WIMP mass with a la statistical uncertainty as ~ 94 GeV ^34% (for a 
WIMP mass of 100 GeV) or ~ 52 GeV+^J (for a WIMP mass of 50 GeV). 

In contrast, since the WIMP scattering spectrum is in principle approximately exponential 
and thus only (very) few WIMP-induced events could be observed in high energy ranges, if we 
have background windows in only high experimental possible energy ranges, the (pretty large) 
contributions from background events could cause (strong) overestimates of the reconstructed 
WIMP masses. It is even worse for large WIMP masses (m x > 100 GeV) 10 . Nevertheless, as 
shown in the lower frame of Figs. 4, with ~ 5% residue background events observed only in 



8 Since for heavy input WIMP masses the reconstructed values are systematically underestimated, probably 
due to the statistical fluctuation with pretty few (~ 50) events discussed later. 

9 Note that here we do not mean that in other energy ranges background events do not exist; in contrast, we 
want to study what could happen once our background discrimination, caused by some natural or even artificial 
reasons, are worse in these energy ranges than others and more background events could thus survive. 

10 Note that the plateau of the lower bound of the statistical uncertainty in the case of a 10% background 
ratio for heavy WIMP masses (m x > 300 GeV) should be caused by our setup for the upper cut-off of the 
reconstructed WIMP mass of 3000 GeV in the simulations. 
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Figure 4: As in Fig. 3, except that the background window in each experiment have been set 
as - 50 keV (upper) and 50 - 100 keV (lower). Note that the background ratios shown here 
are 20% (solid red curves) and 40% (dash-dotted cyan curves) in the upper frame, whereas 5% 
(dotted magenta curves) and 10% (long-dotted blue curves) in the lower frame. 



28 Si + 76 Ge, Q max < 100 keV, Q max bg < 100 keV, exponential bg, 2 x 500 events 
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Figure 5: As in Fig. 3, except that the expected number of total events in each experiment has 
been set as 500. 



high energy ranges, one could in principle still estimate the WIMP mass with a la statistical 
uncertainty as ~ 107 GeV ±Jgj (for an input WIMP mass of 100 GeV) or ~ 58 GeV tf 2 % (for 
an input WIMP mass of 50 GeV). 

Our results shown in Figs. 4 indicate that a small fraction of background events in low energy 
ranges might not affect the reconstructed WIMP mass significantly. However, the WIMP mass 
could be (strongly) overestimated once the same (or even smaller) amount of background events 
exists in high energy ranges. In practice one simple way to reduce the overestimate induced by 
an excess of background events in high energy ranges might be checking the shape of measured 
recoil spectrum. However, considering some suggested modifications of the standard shifted 
Maxwellian velocity distribution, e.g., contributions from discrete "streams" with (nearly) fixed 
velocities [22, 23, 24] or the "late infall" component in the velocity distribution with a velocity 
v ~ v csc [22, 25, 24], it should at least be very careful to reject any recoil event observed in high 
energy ranges artificially. 

In Fig. 5 we rise the expected number of total events in each experiment by a factor of 10, 
to 500 events on average before cuts for the case that residue background events exist in the 
entire experimental possible energy ranges. As shown here, all statistical uncertainties shrink 
by a factor > 3 compared to the results shown in Fig. 3. In addition, the underestimate of 
the reconstructed values of heavy input WIMP masses caused perhaps by the use of pretty few 
(~ 50) events has been reduced with larger data sets; and, the tendency of the underestimate 
of the reconstructed WIMP mass for heavy WIMP masses (m x > 100 GeV) becomes more 
clearly. Finally, Fig. 5 shows that, for the determination of the WIMP mass by using data sets 
of (9(500) total events, the maximal acceptable background ratio could be ~ 10% (i.e., (9(50) 
background events) or even ~ 20%, if WIMPs have a mass of (9(100 GeV). 



3.2.2 Statistical fluctuation 



As discussed in Ref. [6], the statistical fluctuation of the reconstructed WIMP mass by the 
algorithmic procedure in the simulated experiments seems to be pretty problematic, in particular 
for heavier input WIMP masses. Moreover, as mentioned in the previous subsection, with only 
~ 50 total events in each experiment, the tendency of the underestimate of the reconstructed 
WIMP mass for heavier WIMP masses (m x > 100 GeV) seems not to be very clear. Hence, as 
done in Ref. [6], in order to study the statistical fluctuation of the reconstructed WIMP mass 
with different background ratios in our data sets, we consider in this subsection the estimator 
5m introduced in Ref. [6]: 



5m = < 



m xM - m x - m 



if m X) - m < m 
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(27) 
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if m Xjin > m 
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Here m x ; m is the true (input) WIMP mass, m XtTec its reconstructed value, m Xi i i(2) are the 1 (2) a 
lower bounds satisfying x 2 ( m x,i°(i,2)) = X 2 ( m x,rcc) + 1 (4), and m Xi hii(2) are the corresponding 
1 (2) a upper bounds. 

The estimator 5m defined above indicates basically the strength of the deviation of the re- 
constructed WIMP mass from the true (input) value. If the reconstructed la lower and upper 
bounds on the WIMP mass in one simulated experiment cover the true value: m x ^ Q i < m Xjin < 
^x>hii> ^ m * s determined as the deviation of the "reconstructed WIMP mass" from the true 
one in units of the difference between the reconstructed value and the la lower (upper) bound, 
once the reconstructed value is overestimated (underestimated). However, if the true WIMP 
mass lies outside of the experimental la bounds (the reconstructed value is more strongly over- 
/underestimated), 5m is determined as the deviation of the u la lower (upper) bound" from the 
true WIMP mass in units of the difference between the la and 2a lower (upper) bounds. Note 
that, it has been found in Ref. [6] as well as in the results presented in the previous subsection 
that the uncertainty intervals of the median reconstructed WIMP mass are quite asymmetric; 
similarly, the distance between the la and 2a bounds can be quite different from the distance 
between the reconstructed value and the la bound [6]. The definition of 5m in Eq.(27) takes 
these differences into account, and also keeps track of the sign of the deviation: if the recon- 
structed WIMP mass is overestimated (underestimated), 5m is positive (negative). Moreover, 
1 5m | < 1 (2) if and only if the true WIMP mass lies between the experimental 1 (2) a bounds. 

In Figs. 6 we show the normalized distributions of the estimator 5m defined in Eq. (27) for a 
rather heavy input WIMP mass of 200 GeV with 50 (upper) and 500 (lower) total events on av- 
erage before cuts in each experiment. As discussed in Ref. [6], the deviation of the reconstructed 
WIMP mass in the simulated experiments looks asymmetric and non-Gaussian. However, it can 
be seen here clearly that, the more the background events in our analyzed data sets, the more 
concentrated the 5m value in the range between —1 and as well as between and +1. More- 
over, for the case with rather larger data sets of 500 total events, by increasing the background 
ratio the distribution becomes to be more symmetric and Gaussian-like, although the central 
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Figure 6: Normalized distributions of the estimator 5m defined in Eq. (27) for an input WIMP 
mass of 200 GeV with 50 (upper) and 500 (lower) total events on average before cuts in each 
experiment. Parameters and notations are as in Fig. 3. Note that the bins at 5m = ±5 are 
overflow bins, i.e., they also contain all experiments with \5m\ > 5. See the text for further 
details. 



value of 8m seems to fall at ~ —0.5 because of the underestimate of the reconstructed WIMP 
mass. 

In Ref. [6] it has been mentioned that with increasing number of total events the distribution 
of the estimator Sm becomes slowly Gaussian. Figs. 6 here (and Fig. 9 shown later also) indicate 
that with a larger background ratio in the analyzed data sets the distribution of 5m approaches 
to be Gaussian more fast. This interesting observation might be able to offer some new ideas 
for improving the algorithmic procedure for the reconstruction of the WIMP mass with a higher 
statistical certainty. 



3.2.3 With the constant background spectrum 

In order to check the need of a prior knowledge about an (exact) form of the residue back- 
ground spectrum, we consider briefly in this subsection a rather extrem case, i.e., the constant 
background spectrum in Eq. (13). 

In Fig. 7 we show the reconstructed WIMP mass and the lower and upper bounds of the 
lcr statistical uncertainty with mixed data sets as functions of the input WIMP mass. As in 
Figs. 2, the windows of the constant background spectrum are set as the same as the experimental 
possible energy ranges, i.e., between and 100 keV for both experiments. The background ratios 
shown here are no background (dashed green curves), 5% (dotted magenta curves), and 10% 
(long-dotted blue curves) background events in the whole data sets. Each experiment contains 
again 50 total events on average before cuts; all of these events are treated as WIMP signals in 
the analysis. 
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Figure 7: As in Fig. 3, except that the constant background spectrum in Eq. (13) has been 
used. Note that the background ratios shown here are 5% (dotted magenta curves) and 10% 
(long-dotted blue curves). 
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Figure 8: As in Fig. 7, except that the background window in each experiment has been set as 
0-50 keV. 



It can be seen clearly that, as discussed above, since the constant background spectrum has 
relatively flatter shape compared to the WIMP scattering spectrum for not only light, but also 
heavy WIMP masses, and the measured energy spectrum should thus always mimic a scattering 
spectrum induced by heavier WIMPs, the reconstructed WIMP masses are therefore overesti- 
mated for all input WIMP masses, especially for the heavier masses. Actually, the result shown 
here looks more likely that shown in the lower frame of Figs. 4, since in both cases residue 
background events contribute significantly more (compared to the exponential-like WIMP scat- 
tering spectrum) in high energy ranges. Not surprisingly, the larger the background ratio, the 
more strongly overestimated the reconstructed WIMP masses, in particular for the heavier input 
WIMP masses. Nevertheless, for (approximately) constant residue backgrounds with a fraction 
of ~ 5% in background windows as the entire experimental possible ranges, one could in principle 
still estimate the WIMP mass with a la statistical uncertainty as ~ 117 GeV ^35% (for 100 GeV 
WIMPs) or ~ 56 GeV+*jg (for 50 GeV WIMPs). Once WIMPs are light (m x ~ 0(25 GeV)), 
the maximal acceptable background ratio could even be ~ 10%. 

Moreover, as done in Sec. 3.2.1, in Fig. 8 we shrink the background window in each experiment 
to the lower energy range between and 50 keV. Not surprisingly, while for light WIMPs 
(m x < 70 GeV) 11 , relatively more background events still contribute to high energy ranges; 
for heavy WIMPs (m x > 70 GeV), relatively more background events contribute now to low 
energy ranges and, consequently, the reconstructed WIMP masses are therefore underestimated 
for heavy WIMPs. 

11 Remind that the actual value of this "critical" WIMP mass depends in practice strongly on the WIMP 
scattering spectrum as well as on the residue background spectrum and therefore differs from experiment to 
experiment. 




Figure 9: As in the upper frame of Figs. 6, except that the constant background spectrum in 
Eq. (13) has been used. 



'Si + Ge, Q max < 100 keV, Q max bg < 100 keV, constant bg, 2 x 500 events 




10 20 30 50 70 100 200 300 500 700 1000 



m x,in [GeV] 



Figure 10: As in Fig. 7, except that the expected number of total events in both experiments has 
been set as 500. Note that the background ratios shown here are no background (dashed green 
curves), 1% (long-dotted blue curves), and 2% (solid red curves), i.e., a factor of 10 smaller than 
the ratios used before. See the text for further details. 



On the other hand, as in Sec. 3.2.2, in Fig. 9 we check the normalized distributions of the 
estimator 5m for an input WIMP mass of 200 GeV with 50 total events on average before cuts 
in each experiment. It can be seen very clearly that, with increasing background ratio the value 
of 5m concentrates more and more strongly to 2. This means that, due to the contribution from 
residue background events, the reconstructed WIMP mass is most possibly ~ 2a overestimated. 
Moreover, compared to the non-Gaussian form of the distributions for the case with the expo- 
nential background spectrum shown in the upper frame of Figs. 6, the distributions with the 
constant spectrum look more likely Gaussian, despite of the asymmetry due to the overestimate 
of the WIMP mass. Nevertheless, Figs. 6 and Fig. 9 indicate that background events seem to 
let the distribution of the deviation of the reconstructed WIMP mass be more symmetric and 
Gaussian, no matter what kind of energy spectrum they would have. 

Finally, in Fig. 10 we rise the expected number of total events in each experiment by a 
factor of 10, to 500 events on average before cuts for the case that residue background events 
exist in the entire experimental possible energy ranges. Note that the background ratios shown 
here are no background (dashed green curves), 1% (long-dotted blue curves), and 2% (solid red 
curves), i.e., a factor of 10 smaller than the ratios used before. In the lower frame of Figs. 4 
and in Fig. 7, we found that once ~ 5% - 10% events in our analyzed data sets are residue 
backgrounds and (most of) these events are recorded in high energy ranges, no matter what 
kind of spectrum shape they would have, the reconstructed WIMP mass could be (strongly) 
overestimated. However, Fig. 5 and Fig. 10 here show that, by increasing the event number 
and decreasing the background ratio, one could in principle determine the WIMP mass (pretty) 
precisely without knowing the (exact) form of the spectrum of residue background events. 

4 Summary and conclusions 

In this paper we reexamine the model-independent data analysis method introduced in Refs. [5, 
6] for the determination of the mass of Weakly Interacting Massive Particles from data (measured 
recoil energies) of direct Dark Matter detection experiments directly by taking into account a 
fraction of residue background events, which pass all discrimination criteria and then mix with 
other real WIMP-induced events in the analyzed data sets. Differ from the maximum likelihood 
analysis described in Refs. [21, 14, 13], our method requires neither prior knowledge about the 
WIMP scattering spectrum nor about different possible background spectra; the unique needed 
information is the recoil energies recorded in two direct detection experiments with two different 
target nuclei. 

In Sec. 2 we considered first the measured energy spectrum for different WIMP masses with 
two forms of possible residue background spectrum: the target-dependent exponential spectrum 
and the constant spectrum. The exponential background spectrum contributes relatively more 
events to high energy ranges once WIMPs are light (m x < 100 GeV), and to low energy ranges for 
heavy WIMP masses (m x > 100 GeV); whereas the constant background spectrum contributes 
always relatively more events to high energy ranges. As the consequence, the energy spectrum of 
all observed events looks more likely to be a scattering spectrum induced by heavier WIMPs, once 
the spectrum of residue background events (induced perhaps by two or more different sources) is 
either exponential-like (and WIMPs are light) or approximately constant (for all WIMP masses); 
while if WIMPs are heavy and the residue background spectrum is approximately exponential, 
the measured energy spectrum would look more likely to be a scattering spectrum induced by 
lighter WIMPs. 

In Sec. 3.2 the data sets generated in Sec. 2 have been analyzed for reconstructing the mass 



of incident WIMPs by using the model-independent method. With the exponential background 
spectrum, the input WIMP mass would be overestimated once WIMPs are light (m x < 100 GeV), 
or, in contrast, would be underestimated for heavy WIMPs (m x > 100 GeV). Our simulations 
show that, for background windows in the entire or low experimental possible energy ranges, 
one could in principle reconstruct the WIMP mass with a maximal fraction of ~ 20% of residue 
background events in the analyzed data sets; whereas for background windows in high energy 
ranges, the maximal acceptable fraction of residue backgrounds is only ~ 10%. 

Moreover, in order to check the need of a prior knowledge about an (exact) form of the 
residue background spectrum, we considered also the case with the constant background spec- 
trum. In this rather extrem case, the WIMP mass would always be overestimated, especially 
for heavy WIMPs (m x > 100 GeV). Our simulations give then a maximal acceptable fraction 
of ~ 5% - 10% of residue background events in the data sets for background windows in the 
entire or low experimental possible energy ranges. Nevertheless, we found also that, by means of 
increased number of observed (WIMP-induced) events and improved background discrimination 
techniques [9, 11], the WIMP mass could in principle be determined (pretty) precisely, no matter 
what kind of energy spectrum residue background events would have. 

On the other hand, in order to check the statistical fluctuation of the reconstructed WIMP 
mass with increased background ratio, we considered also the distribution of the deviation of 
the reconstructed WIMP mass from the true value. It was found in Ref. [6] that, for a rather 
heavy WIMP mass of 200 GeV, the distribution of the deviation of the reconstructed WIMP 
mass is asymmetric and non-Gaussian, either with data sets of only a few (0(50)) events or with 
larger date sets (of O(500) events). However, our simulations with different background ratios 
show that, firstly, for both used (exponential and constant) background spectra, with increasing 
background ratio the distribution of the deviation of the reconstructed WIMP mass becomes 
more and more concentrated, although still asymmetric and non-Gaussian. Secondly, for the 
more realistic exponential background spectrum and using data sets with a larger number of total 
events, with increasing background ratio the distribution of the deviation becomes somehow more 
symmetric and Gaussian. This observation might be able to offer some new ideas for improving 
the algorithmic procedure for the reconstruction of the WIMP mass with a higher statistical 
certainty. 

In summary, our study of the effects of residue background events in direct Dark Matter de- 
tection experiments on the determination of the WIMP mass shows that, with currently running 
and projected experiments using detectors with 10 -9 to 10 -11 pb sensitivities [26, 27, 10, 28] and 
< 10~ 6 background rejection ability [9, 11, 12, 8], once two or more experiments with different 
target nuclei could accumulate a few tens events (in one experiment), we could in principle al- 
ready estimate the mass of Dark Matter particle with a reasonable precision, even though there 
might be some background events mixed in our data sets for the analysis 12 . Moreover, two forms 
for background spectrum and three windows for residue background events considered in this 
work are rather naive. Nevertheless, one should be able to extend our observations/discussions 
to predict the effects of possible background events in their own experiment. Hopefully, this will 
encourage our experimental colleagues to present their (future) results not only in form of the 
"exclusion limit(s)", but also of the "most possible area(s)" on the cross section versus mass 
plan. 

12 A possible first test could be a combination of the events observed by the CoGeNT experiment with their 
Ge detector with the events observed in the oxygen band of the CRESST-II experiment [29, 30, 31]. 
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A Formulae needed in Sec. 3 

Here we list all formulae needed for the model-independent data analyses procedure used in 
Sec. 3. Detailed derivations and discussions can be found in Refs. [7, 6]. 

A.l Estimating r(Q m - in ), I n (Qmim Qmax)? and their statistical errors 

Firstly, consider experimental data described by 

Qn- b T<Qn,i<Qn+ b T, * = 1, 2, • • • , N n , 71 = 1, 2, • • • , B. (Al) 

Here the entire experimental possible energy range between Qmin and Q max has been divided 
into B bins with central points Q n and widths b n . In each bin, N n events will be recorded. Since 
the recoil spectrum dR/dQ is expected to be approximately exponential, the following ansatz 
for the measured recoil spectrum (before normalized by the exposure £) in the nth bin has been 
introduced [7]: 

fdR\ =fdR\ = rn e fe ^-Q s ,n). (A2 ) 

V^/expt, n \dQJ expt, Q~Q n 

Here r n is the standard estimator for (dR/dQ) cxpt at Q = Q n : 
N 

r» = > (A3) 

k n is the logarithmic slope of the recoil spectrum in the nth Q— bin, which can be computed 
numerically from the average value of the measured recoil energies in this bin: 

Q^Q- n \ n = (|) coth - i- , (A4) 

where 

i N„ 

(Q - Qn) X \n = W E (Qn,i ~ Qn)" ■ (A5) 
iV « i=l 

The error on the logarithmic slope k n can be estimated from Eq. (A4) directly as 



Akn) =kUl 



k n b n /2 



sinh(/c n 6 n /2) 



2 



-2 



(? 2 (Q-Qn\n), (A6) 



with 



a 2 (Q - Q n \n) 



N„ 



(Q Qn) 2 \n Q Qn\ n 



(A7) 



Q St n in the ansatz (A2) is the shifted point at which the leading systematic error due to the 
ansatz is minimal [7], 



Qs,n Qn 4" 7 hi 



sinh(fc ra 6 n /2) 
k n b n /2 



(A8) 



Note that Q s , n differs from the central point of the nth bin, Q n . From the ansatz (A2), the 
counting rate at Q — Q m m can be calculated by 



r(Q mi n) = ne* 1 ^"^ 1 ' , 
and its statistical error can be expressed as 



cr 2 (r(Q miQ )) = r 2 (Q min ) 



1 



+ 



1_/V 

h 1 2 



1 + coth 



<J 2 (h) 



since 



2/ \ N n 

° (r n ) = -rr ■ 

Finally, since all I n are determined from the same data, they are correlated with 

Q(n+m~2)/2 



cov(J n , I m ) = 



t F\Qa) ' 



(A9) 



(A10) 



(All) 



(A12) 



where the sum runs again over all events with recoil energy between Qmin and Qmax- And the 
correlation between the errors on r(Q min ), which is calculated entirely from the events in the 
first bin, and on I n is given by 



cov(r(Q min ),I n ) 

f(Qmin) In(.Qmim Qm'm + 6i) 



x < — + 



x 



1 + coth 



h v 2 J V ' V 2 

In+2(Qmm,Qmm + h) _ J_ _ f bj\ ( b X k X 

L /n(Qmin,Qmin + &l) ^ h \ 2 ) °° ^2, 



a 2 (h) ; (A13) 



note that the sums Jj here only count in the first bin, which ends at Q — Q m i n + b\. 



A. 2 Statistical errors on m x given in Eqs. (18) and (21) 

By using the standard Gaussian error propagation, a lengthy expression for the statistical error 
on m x \^ given in Eq. (18) can be obtained as 



XJ\(v n ) 



^m x /m Y \m x - m Y \ (K^x/Kny) 



(jln,x/'R-n,Y - y/mx/my) 



X 



E 



'dU n . 



X 



'&JZ, 



n,X 



dc 



cov(q iX , c j>x ) + (A — >Y) 



1/2 



• (A14) 



Here a short-hand notation for the six quantities on which the estimate of m x depends has been 
introduced: 



Cl,X — I n ,x , 



C2,X — h,x , 



C3,X — r x(Qmm,x) 



(A15) 



and similarly for the c^y. Estimators for cov(q,c j ) have been given in Eqs. (A12) and (A13). 
Explicit expressions for the derivatives of 1Z UyX with respect to c^ x are: 



dU 



n,X 



dl n , x 
dTZ ni x 



n+ 1 

n 

1 

n 



F x (Qmin,x) 



2gS /2 ^(g m in,x) + (n + i)/n,xFKg min ,x) 



X 



F X {Q 



min,X, 



2 <5L / in,X r x(Qmin,x) + h,X Fx(Qmm,x) 



K 



n,X 



(Al6a) 



(Al6b) 



and 



dTZ 



n,X 



dr x (Q 



min,X J 



2 

n 



n (n+l)/2 j 



(n + l)Q 



1/2 . 
min,X I n,X 



2Q ( Xx 2 rx(Qnun,x) + + 1)/„,xF£ (Q m in,x) . 

min,x) 



X 



F X (Qr 



1 /2 

2 <5min,X r x(Qmin,x) + h,xF x (Q m m,x) 



K 



n,X 



(A16c) 



explicit expressions for the derivatives of 1Z n y with respect to Qy can be given analogously. Note 
that, firstly, factors TZ n ,(x,Y) appear in all these expressions, which can practically be cancelled 
by the prefactors in the bracket in Eq. (A14). Secondly, all the Iq<x,y) and I n <x,Y) should be 
understood to be computed according to Eq. (17) with integration limits Q m i n and Qmax specific 
for that target. 

Similar to the analogy between Eqs. (18) and (21), the statistical error on m x \ a given in 
Eq. (21) can be expressed as 



KK) 5/2 


m x — my 


(K a , x /K a , Y ) 




Tla,x/Tla,Y - (m x /m Y ) 5/2 


2 



X 



1 3 
IT 2 ^ 

K-ajc i,j=2 



dn 



a,X 



X 



da 



i,X 



dcj 



cov(c iiX , c jtX ) + (X — >Y) 



1/2 



, (A17) 



where we have again used the short-hand notation in Eq. (A15); note that Ci^ x ,y) = ^n,(x,Y) 
does not appear here. Expressions for the derivatives of 1Z a ,x can be computed from Eq. (22) as 



dTZ a . 



x 



dl ,x 

on a>x 

dr x (Q 



F 2 X (Q 



min,Jf / 



2Qmin,X r x(Qmm,x) + Iq,xF x 

(Qmin,x) 



a,X , 



min,X, 



2Q 



1/2 
min,X 



2Qmin,X r x{Qmm,x) + h,xF x (Qram,x) 



K 



<T,X 



(Al8a) 



(Al8b) 



and similarly for the derivatives of lZ a y. 



A. 3 Covariance of fa defined in Eqs. (24a) and (24b) 

The entries of the C matrix given in Eq. (25) involving basically only the moments of the WIMP 
velocity distribution can be read off Eq. (82) of Ref. [7], with an slight modification due to the 
normalization factor in Eq. (24a) 13 : 



cov (fa fj) = M x 



Here we have defined 



a 



and 



Di = 



300 km/s ' 



U f 3 cov(/ , J ) + + + l)cov(4 /,-) 

- &{j + cov(/ , Ij) - a\i + l)fj cov(J , h) 

+ D i D j a 2 (r(Q niin )) - (DJj + Djfc) cov(r(Q min ), 7 ) 

+ a j (j + 1) A cov(r(g min ), Ij) + a\i + l)Dj cov(r(Q min ), h) 

(A19) 



(A20) 
(A2I) 



'min^ '{Q min) / ' F^iQr, 

a 



dr(Q min ) 



" Vmin Vmin H I i 



for « = —1, 1, 2, . . . , n max ; and 
2 



n = f-n 1/2 f \ 

^"max+l F 2 (Q ) ^ ^minjnmax+iy ■ 



(A22a) 



(A22b) 
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